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We solve the elliptic equations associated with the Hamil- 
tonian and momentum constraints, corresponding to a system 
composed of two black holes with arbitrary linear and angu- 
lar momentum. These new solutions are based on a Kerr- 
Schild spacetime slicing which provides more physically re- 
alistic solutions than the initial data based on conformally 
flat metric/maximal slicing methods. The singularity/inner 
boundary problems are circumvented by a new technique that 
allows the use of an elliptic solver on a Cartesian grid where 
no points are excised, simplifying enormously the numerical 
problem. 

PACS number(s): 04.70.Bw,04.25.Dm 



I. INTRODUCTION AND METHOD 

Any code designed to evolve a general relativistic sys- 
tem will start with a solution to the initial value problem 
corresponding to the astrophysical situation of interest. 
In the 3+1 formulation of general relativity, the problem 
is defined by the Hamiltonian and momentum constraints 






(1) 



Above, R is the 3-dimensional Ricci scalar constructed 
from the spatial 3- metric gij, and K and Aij are the 
trace and trace-free parts of Kij such as 
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The problem of solving the Hamiltonian and momentum 
constraints for two black holes has been addressed in the 
past by several groups (see [0 and references therein). 
For us, it is inherently 3-dimensional since we want to 
be able to specify arbitrary boost and spin directions for 
each hole. The methods in most frequent use MM are 
based on an approach which chooses maximal spatial hy- 
persurfaces {K = 0) and takes the spatial 3-metric to be 
conformally flat {gij = (fi'^Sij, being the conformal fac- 
tor). Under these conditions, the Hamiltonian constraint 
can be decoupled from the momentum constraint. Ana- 
lytical solutions for the extrinsic curvature Kij for holes 
with specific linear momenta can be found 0] . This sim- 
plification leaves only one elliptic equation for 0, which 



is derived from the Hamiltonian constraint. The inner 
boundaries (the throats of the black holes) are usually 
dealt with by imposing an isometry condition between 
two identical asymptotically flat spatial slices, joined by 
an Einstein-Rosen bridge, though other boundary condi- 
tions are sometimes used. Unfortunately, the numerical 
solution of the equation for presents a technical chal- 
lenge at the inner boundaries. Brandt and Briigmann 
Ig simplified this problem by compactifying the internal 
asymptotically flat regions to obtain a domain without 
inner boundaries. However, the main disadvantage of 
these methods is not numerical, but related to the physi- 
cal interpretation of black-hole spaces described through 
a conformally flat 3-metric. There are no space slices 
for which the spatial 3-metric of a single Kerr (non-zero 
spin) black hole can be written in a conformally flat way 
M , and recent work on sequences of initial data sets for 
circular orbits casts some doubt on the physical realism 
of conformally flat approaches to black hole binaries [Q . 
One way to overcome this problem is to use a Kerr-Schild 
slicing of spacetime [g. Matzner, Huq, and Shoemaker 
Q proposed a method that bases the initial data on a 
background metric and extrinsic curvature that is a su- 
perposition of Kerr black hole metrics written in ingo- 
ing Eddington-Finkelstein coordinates. Thus, for a sin- 
gle Kerr black hole one obtains the exact solution to the 
problem. Marronetti et al. 10 1 added to this method a 



variation on the background flelds that eliminates the in- 
ner boundary problem, greatly simplifying the numerical 
treatment of the elliptic equations. The background ^e\As 
generated in this way are also very good approximate so- 
lutions to the initial value problem (i.e., the violation of 
the constraints (|l|) is small enough to fall below the nu- 
merical truncation error for a wide range grid spacings). 

We present here the first solution to the initial value 
problem for black-hole binaries based on physically re- 
alistic background 3-metric and extrinsic curvature. We 
study the particular case of a hyperbolic encounter of two 
Kerr black holes of mass m, separated by a distance of 
11.5 m. 

The method begins by specifying a conformal spatial 
metric which is a straightforward superposition of two 
Kerr-Schild single hole (spatial) metrics: 

fjij = Sij + 2 iB iH ili ilj + 2 2R 2H 2U 2ij , 

K^iB iK^'+2B2K,' , 

1 



A., =5n(. iiB,K^^ + 2B2K,;~^S^;K) , (2) 

where the parenthesis in the subscripts denote sym- 



metrization in i and j. The scalar function H has a 
known form B , and Ix is an ingoing null vector congru- 
ence associated with the solution. The fields marked with 
the pre-index 1 (2) correspond to an isolated black hole 
with specific angular momentum ai (32) and boosted 
with velocity vi (V2). The trace- free part of the ex- 
trinsic curvature tensor Aij is also constructed from the 
superposition of the curvatures iK^ and 2K^ associated 
with the black holes 1 and 2 respectively. The atten- 
uation function iB {2B) is unity everywhere except in 
the vicinity of hole 2 (hole 1) where it rapidly vanishes, 
so the fields there are effectively those of a single black 
hole, thus providing an exact solution to the constraints 
for distances arbitrarily close to the singularities llfl] . 

Following the conformal decomposition presented by 
York and collaborators |jll| , we relate the physical metric 
gij and the trace-free part of the extrinsic curvature Aij 
to the background fields through a conformal factor: 



^4- 



9t3 = </> 

A'^ = 6-^° A'' 



(3) 



In order to find a solution to the four constraint Eqs. (|^), 
we add a longitudinal part to the extrinsic curvature A''^ : 



A*^=0-i°(i*^' + (?u;)') , (4) 

where w* is a vector potential to be solved for and 

2 



(Iw) 



V'w^ 






(5) 



Plugging Eqs. (|3;5|), into the Hamiltonian and momen- 
tum constraints ([1 ) , we obtain four coupled elliptic equa- 
tions for the fields and ui* ||ll| : 



r\A''^ + {iwr){A, + {iw),,)) 






(6) 



To solve the elliptic Eqs. (O) we use an adaptation 
of an multigrid elliptic solver developed for the solution 
of the initial value problem of neutron-star binaries p2[ , 
which uses second order stencils. To optimize the use of 
memory, this solver was designed to handle only identical 
black holes; only one object is placed on the numerical 
grid, the second is accounted for using reflective bound- 
ary conditions at the corresponding surface of the grid 
cube (see 111] for details) . The four elliptic equations are 
solved consecutively in an iteration loop. The iteration 
starts with a trivial initial guess, namely 



</>(x*) = 1 
w\x')^ {0,0,0) \/x' 



(7) 



and it relaxes the fields until the variation from one it- 
eration step to the next falls below some fraction of the 
truncation error. This typically takes around 10 to 20 it- 
eration steps. The sources for these elliptic Eqs. consist 
mostly [Of of the residuals of the constraint equations as 
defined in llQ]. The use of attenuation functions elimi- 
nates the singular behavior of these sources near the ring 
singularities, simplifying the numerical implementation 
of the elliptic solver. Instead of using excision techniques 
around the singularities [M, we modified a standard el- 
liptic solver to "ignore" (i.e. leave the fields with the 
initial values (^) the grid points in a small volume em- 
bedding the singularity. This "inner" region is defined 
monitoring the value of the background 3-dimensional 
Ricci scalar |1J], and is chosen to be smaller than the 
excision volume used in current numerical simulations 
[]l4| , making the initial value data suitable for the evo- 
lutionary codes. The choice of the "inner" region values 
(0) is based on the fact that arbitrarily close to the ring 
singularity, the background fields recover the single-hole 
values which are an exact solution to Eqs. (P). We use 
Robin conditions for the outer boundaries jO], guaran- 
teeing that 4) and w^ fall off as constant + 0{\/r) at 
infinity. 

The calculations presented in this paper were per- 
formed using 16 processors of an Origin2000 system lo- 
cated at the National Computational Science Alliance 
center in Urbana- Champaign, taking close to 200 CPU 
hours. 



II. RESULTS AND CONCLUSIONS 

The method described above was applied to the 
problem of two black holes with mass m in a hyper- 
bolic encounter configuration with parameters ri — 
(5.75 m,0,0), vi = (0,0.5 c,0), and ai = (0,0,0.5 m) 
corresponding to black hole 1 coordinate position, veloc- 
ity, and specific angular momentum. The parameters for 
black hole two are y-2 = — fi, V2 — — Vi, and a2 = ai. 
To test the convergence of the code, two grid resolutions 
were used: one with grid spacing h — to/4 (low) and an- 
other with h = to/8 (high). The number of grid points 
was 60,750 and 536,238 respectively and the bound- 
ing box of the numerical grid covers a rectangle from 
(—11.5 TO, —5.75 771, —5.75 to) to (11.5 to, 5.75 to, 5.75 to,). 
We estimated the ADM mass corresponding to the back- 
ground fields and to the full numerical solution, obtaining 



2.34 TO (note that (1 



-1/2 = 1.155) and 2.18 to re- 



spectively. This was done by integrating over the bound- 
ary surface of the computational grid. The difference 
is accounted for by the binding energy. We also ob- 
tained a background total angular momentum estimation 
of 7.99 nn? , consistent with analytical estimates. 

When differenced under our numerical scheme, even 
analytically exact single-hole solutions show inevitable 



nonzero residual errors in the Hamiltonian and momen- 
tum constraints, which increase near the hole, due to 
the larger gradients there. In Figures || and g these are 
plotted as points (filled circles for resolution to/8, filled 
squares for resolution to/4). The decrease of this trunca- 
tion error with increasing resolution shows the convergent 
behavior of the differencing scheme. Also in Figures 1-2 
we plot as curves the computed constraint residuals of the 
rightmost hole in our two-hole solution. Except for the 
innermost points in the solution, the residuals are effec- 
tively identical to those of the analytical single-hole solu- 
tion at the same resolution. The inset in Figure El shows 
this discrepancy at the innermost points for the high reso- 
lution curve (empty squares mark the grid points and the 
singularity is represented by a gray vertical band). We 
thus confirm our two-hole solution at the level of trunca- 
tion error except for these innermost points, which will 
anyway be excised in the evolution scheme n& . 



the great degree of accuracy of the background fields as 
approximate solutions to the constraints pO| . 
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FIG. 1. Residuals of the Hamiltonian constraint near the 
hole located at (5.75 m,0,0). The full (dashed) line corre- 
sponds to the two-hole calculation done with grid spacing 
h = m/8 {h = m/4). The circles (squares) correspond to the 
computation of the Hamiltonian constraint of a single black 
hole on a grid with spacing h = m/8 (h = m/4). 

Figure ^ is a contour plot of the conformal factor (j). 
The inner thick lines limit the "inner" region where the 
initial value (j) — 1 is preserved. The apparent horizons 
of the holes were calculated using a code developed by 
Huq et al. |17|| and are represented by the outer thick 
lines. The resulting horizon areas for the background 
system and the full solution were 2 x 39.9 m? and 2 x 
38.4 m? respectively, which are consistent with the value 
of 45.7 m? obtained for a single- hole configuration . 

Table | shows the li (average of the absolute value) 
and loo (maximum value) norms for the calculated fields. 
The final fields (j) and w' depart very little from the ini- 
tial guess (R) all across the numerical grid, indicating 



FIG. 2. Residuals of the Momentum constraint (x comp.) 
near the hole located at (5.75 m,0,0). The notation is the 
same as in Figure hi 

The attenuation method, while providing a clear nu- 
merical advantage, also raises the question of physical in- 
terpretation of these results. One would expect the fields 
around each hole to be tidally distorted by the presence 
of the second hole, and the attenuation functions par- 
tially occlude this effect. However, we notice that most 
of this effect is confined to the regions very near the ring 
singularities: i.e. regions that will be excised from the 
computational domain by the time evolution code. At 
the same time, we see from figure || that the tidal inter- 
action is still (at least partially) present at the location of 
the apparent horizons. The influence of the attenuation 
functions can be interpreted as an addition of gravita- 
tional radiation that conveniently (from the numerical 
point of view) eliminates the singular behavior from the 
elliptic equations near the singularities. 

This solution to the initial value problem presents a dif- 
ferent approach than the conformally flat/maximal slic- 
ing methods, and allows us to specify more directly the 
physical content of the data. It also entails a compu- 
tational simplification on the treatment of the singulari- 
ties on the working grid. However, the physical aspects 
of these data sets can be further refined, for instance 
by providing background fields based on Post-Newtonian 
expansions or with a more physical prescription for the 
background fields for circular orbits [Q . 
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FIG. 3. Contour plot of the conformal factor </>. The inner 
thick lines limit the "inner" region where the initial value 
(j) = 1 is preserved. The outer thick lines show the apparent 
horizons and the arrows the direction of the boosts. 
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TABLE I. Maximum (loo) and average of the absolute val- 
ues (li) of the conformal factor and vector potential w^ on 
the grid volume. 



l.OOOe 4- 
0.997e + 



0.462e 
0.745e 



0.168e 
0.146e 



0.421e - 2 
0.452e - 3 



